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Abstract 


Wc  present  new  results  from  a  joint  ocean -acoustic  modeling  study  of  solitary  wave  generation  in  the  Strait  of  Messina, 
their  propagation  in  the  Tyrrhenian  Sea  and  subsequent  shoahng  in  the  Gulf  of  Gioia.  The  non  hydro  static  3D  EULAG 
model  is  used  for  the  oceanographic  predictions.  The  simulations  are  initialized  with  measured  temperature  and  salinity 
profiles  from  an  October  1995  survey  of  the  Messina  region,  and  forced  with  semidiurnal  tidal  magnitudes  predicted  by 
a  barotropic  tidal  model.  Parameter  sensitivity  studies  are  performed.  The  predicted  solitary  wave  trains  are  compared 
with  CTD  chain  measurements.  The  model  results  and  data  are  examined  through  a  wavelet  analysis.  The  wavelengths 
are  tracked  by  the  spines  {maximum  intensity  for  each  wavelength)  at  various  times.  From  the  slope  of  the  variations, 
phase  speeds  are  derived  as  a  function  of  wavelength.  For  the  parameters  extracted  from  CTD  measurements  and  existing 
tidal  conditions,  phase  speed  distribution  for  wavelengths  ranging  from  about  0.6  m  to  1 .6  km  are  obtained.  The  model 
predicted  phase  speed  magnitudes  range  from  0,85  m  s_1  to  0.93  ms"1.  The  phase  speeds  derived  from  data  range  from 
0,77  m  s  1  to  0.88  m  s  1 .  The  model  predicted  phase  speed  versus  wavelength  distribution  has  similar  trends  to  the  phase 
speed  versus  wavelength  distribution  derived  from  data.  The  shoaling  of  the  solitary  waves  in  the  Gulf  of  Gioia  is  studied. 
Calculations  of  the  acoustical  field  are  conducted,  along  the  solitary  wave  propagation  path,  with  the  parabolic  (PE)  acous¬ 
tical  model. 

Published  by  Elsevier  Ltd. 

Keywords:  Ocean  models;  Solitary  waves;  Barotropic  tide;  Messina  sill;  Wavelet  analysis;  CTD  chain  data 


L  Introduction 

The  Strait  of  Messina  connects  the  Tyrrhenian  Sea,  on  the  north  side,  and  the  Ionian  Sea,  on  the  south  side. 
Fig.  1  shows  the  area  of  interest,  with  bathymetric  contours  in  meters.  The  strait  contains  a  sill  that  raises  to 
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Fig,  I  Strait  of  Messina  and  Gulf  of  Gioia  region.  Triangles  are  CTD  stations.  CTD  chain  tracks  5,  6,  and  8  are  indicated  (green  curve). 
The  dashed  box  shows  the  horizontal  model  domain. 


within  80  m  of  the  surface  at  the  shallowest  point.  The  water  masses  of  the  Ionian  and  Tyrrhenian  Seas  mix  in 
the  strait.  On  the  north  side  of  the  strait  (towards  the  Gulf  of  Gioia),  the  depth  increases,  reaches  a  maximum 
and  then  decreases. 

One  of  the  first  attempts  to  explain  the  tidal  dynamics  of  the  Straits  of  Messina  was  by  Sterneck  (1915).  He 
proposed  a  pattern  of  tidal  oscillations,  in  the  Strait  of  Messina,  generated  by  the  M2  semi-diurnal  tide  and 
obtained  numerical  solutions  of  the  linear  tidal  equations.  Vcrcelli  (1925)  conducted  a  survey  of  the  Strait  ot 
Messina  and  concluded  that  the  current  distribution  is  largely  due  to  tidal  action,  He  showed  that  there  is  a 
large  gradient  of  tidal  amplitude  between  the  tides  north  and  south  of  the  sill  because  the  tides  are  out  of  phase 
by  about  180°.  Because  of  the  phase  opposition  and  the  topographic  constrains  the  tidal  currents  can  reach 
magnitudes  of  around  3ms  '.  The  magnitudes  of  these  currents  were  know  in  ancient  times,  dating  back 
to  Homer’s  Odyssey  (Alpcrs  and  Salusli,  1983).  The  advent  of  the  SEASAT  satellite  has  brought  forth  new 
information  on  the  oceanography  of  the  strait.  Analysis  of  the  SEASAT  SAR  data  showed  the  presence  of 
internal  waves  that  were  linked  to  tidal  currents  moving  over  the  Messina  sill  (Alpers  and  Salusti,  1983). 

The  generation  and  subsequent  propagation  of  solitary  waves  over  topographic  changes  -  c.g.,  at  shelf- 
breaks  or  sills  -  has  been  modeled  with  nonhydrostatic  hydrodynamic  models.  The  hydrodynamic  modeling 
of  tidal  flow  over  steep  topography  shows  that  the  tidal  flow  depresses  the  pycnocline  and  generates  an  inter¬ 
nal  bore  (Lamb.  1994;  Brand  et  al.,  1997;  Warn-Vamas  et  al.,  2003).  The  internal  bore  propagates  and  its  lead¬ 
ing  edge  steepens  through  nonlinear  effects.  Then  frequency  and  amplitude  dispersion  sets  in,  and  the  leading 
edge  disintegrates  into  propagating  solitary  waves  (Lamb,  1994;  Brand  et  al,,  1997;  Warn-Varnas  et  al.,  2003). 
in  particular,  the  interaction  of  the  Strait  of  Messina  sill  with  the  semidiurnal  tidal  motion  results  in  the 
formation  of  internal  bores.  As  the  semidiurnal  tidal  motion  reverses  itself,  the  internal  bores  undergo  a 
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hydraulic  jump  over  the  sill.  Internal  bore  depressions  are  most  pronounced  in  the  thermocline  and  halocline 
regions.  During  the  jump  of  an  internal  bore  over  the  sill,  temperature  and  salinity  fronts  are  formed.  These 
fronts  separate  the  water  masses  in  the  lower  and  upper  parts.  The  interface  between  the  two  water  masses  has 
been  investigated  in  Hopkins  et  al  (1984)  and  Del  Riceo  (1981). 

The  propagating  solitary  wave  trains  are  seen  as  evidence  of  an  isopycnal  depression  (Osborne  and  Burch. 
1980).  The  depressions  associated  surface  convergence  and  divergence  patterns  have  been  observed  in  SAR 
imagery  (Alpcrs  and  Salusti,  1983).  The  isotherm  displacements  can  be  measured  with  a  towed  CTD  therm¬ 
istor  chain  (Sellsehopp,  1997).  The  propagating  solitary  wave  trains  undergo  dispersion  in  amplitude,  wave¬ 
length,  and  phase  speed  (Brand  ct  al.,  1997;  Wam-Varnas  et  al.,  2003).  On  the  north  side  of  the  sill  the  internal 
solitary  waves  propagate  towards  the  Gulf  of  Gioia.  In  the  Gulf  of  Gioia  the  solitary  waves  shoal  as  the  topog¬ 
raphy  raises  toward  the  surface.  The  JANE  1984  (Guardian!  et  al.,  1988;  JANE,  1985)  cruise  measurements 
indicated  a  patch  of  mixed  water  between  40  m  and  100  m  that  could  be  caused  by  mixing  due  to  the  shoaling 
of  internal  solitary  waves. 

In  this  study,  we  model  numerically  the  three-dimensional  oceanographic  structure  from  the  Messina  Strait 
sill  to  the  Gulf  of  Gioia.  This  encompasses  the  channel  for  the  Messina  sill  into  the  Ionian  Sea  and  the 
Tyrrhenian  Sea  into  the  Gulf  of  Gioia.  Our  objectives  are  to  study  the  characteristics  and  effects  of  solitary 
waves  embedded  in  a  three-dimensional  structure  during  generation,  propagation,  and  shoaling.  The  numer¬ 
ical  model  is  initialized  from  a  survey  of  the  region,  conducted  in  October  1995,  and  the  characteristics  of  the 
predicted  solitary  trains  are  compared  against  CTD  chain  measurements.  Effects  of  oceanographic  parameter 
variations  on  solitary  wave  characteristics  are  studied.  The  variations  of  oceanographic  parameters  are  related 
to  the  characteristics  of  solitary  wave  trains  measured  with  the  CTD  chain.  The  shoaling  of  the  solitary  waves 
in  the  Gulf  of  Gioia  is  studied. 

The  temporal  and  spatial  scales  of  tidal ly  induced  internal  bores  and  solitary  waves  are  such  that  they  can 
have  a  significant  effect  on  the  acoustic  field  through  the  sound  speed  structure.  At  certain  frequencies  the 
interaction  of  the  acoustic  field  with  the  solitary  waves  can  be  quite  significant.  Many  interactions  of  the 
acoustical  field  with  the  solitary  wrave  train  can  occur  between  the  source  and  the  receiver  (Warn-Varnas 
et  al.,  2003). 

The  paper  is  organized  as  follows.  The  following  section  briefly  characterizes  the  data  that  drive  the  numer¬ 
ical  model  summarized  in  Section  3  and  Appendix.  Section  4  discusses  ocean  dynamics  simulated  by  the 
model;  whereas  the  validation  of  the  model  results,  using  comparisons  with  measurements  and  parameter- 
sensitivity  study,  is  presented  in  Section  5.  Section  6  discusses  acoustical  effects,  and  the  summary  of  our 
conclusions  is  presented  in  Section  7. 


Temperature  (*C)  Sal  i  nily  (ppt) 

Fig,  2,  Measured  (blue)  and  lined  (red  and  green;  case  I  and  3  in  Table  1  of  Section  4.1.  respectively)  ambient  profiles  of  temperature  (left 
plate)  and  salinity  fright  plate). 
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2 .  Data 

In  October  of  1995  a  survey  of  the  Atlantic  Ionian  Stream  was  conducted  in  the  Strait  of  Sicily*  At  the  end 
of  this  sea  trial,  on  October  24  and  25,  solitary  wave  trains  were  tracked  in  the  Strait  of  Messina  with  a  towed 
conductivity- temperature-depth  (CTD)  chain  (Sellschopp,  1997),  The  chain  had  83  sensors  attached  on  a 
270  m  cable.  The  accuracies  of  the  sensors  were  0*01  K  for  temperature  and  0,02  ppt  for  salinity.  The  resultant 
horizontal  and  vertical  resolution  was  5  m  and  2.5  m  at  a  ship  speed  of  2.5  m  s  1  (Sellschopp,  1997). 

On  the  Tyrrhenian  Sea  side,  tows  5,  6,  and  8  (shown  as  green  lines  in  Fig.  1)  passed  through  the  same 
packet.  The  CTD  chain  measurements  were  corrected  for  ship  motion.  From  the  measured  temperature 
and  salinity  values.  Fig.  2,  the  densities  and  sound  speeds  were  calculated. 

3.  Model 

3.1.  Overview 

The  numerical  predictions  were  conducted  with  the  nonoscillatory  forward-in-time  (NFT)3  model,  termed 
EULAG  for  its  capability  to  solve  the  fluid  equations  in  either  an  Eulerian  (flux-form)  or  a  Lagrangian  (advec- 
tive  form)  mode  (see  Smolarkiewicz  and  Prusa,  2002,  for  a  succinct  review').  The  default  analytic  formulation 
of  EULAG  assumes  the  non  hydrostatic  anelastic  equations  of  motion,  with  options  available  for  comprcss- 
ible/incompressible  Boussinesq,  and  incompressible  Euler/Navier-Stokes  equations  (Smolarkiewicz  et  al, 
2001  f  The  ocean  derivative  of  EULAG  developed  for  the  purpose  of  this  study  builds  on  the  classical  incom¬ 
pressible  Boussinesq  approximation  that  enables  efficient  semi-implicit  integrations  of  stiff  governing  equa¬ 
tions  supporting  gravity  weaves  on  a  broad  range  of  scales.  The  resulting  NFT  ocean  model  retains  all 
multiscale  benefits  of  the  EU  LAG’s  mathematical/numerical  design,  widely  documented  in  the  literature. 
Among  others,  these  include:  (i)  formulation  and  solution  of  governing  PDEs  in  generalized  time-dependent 
curvilinear  coordinates,  admitting  grid  adaptivity  to  flow'  features  and/or  irregular  boundaries  (Prusa  and 
Smolarkiewicz,  2003;  Wedi  and  Smolarkiewicz,  2004;  Smolarkiewicz  and  Prusa,  2005;  Prusa  and  Gutowski, 
2006);  and  (ii)  the  direct  numerical  simulation  (DNS),  large-eddy  simulation  (LES),  and  implicit  large-eddy 
simulation  ( I  LES)  turbulence  modeling  capabilities,  facilitating  applications  al  broad  range  of  Reynolds  num¬ 
bers  (Smolarkiewicz  and  Prusa,  2002;  Smolarkiewicz  and  Margolin,  2007). 

For  clarity,  we  present  only  the  adiabatic  meso-scale  model  equations,  and  dismiss  the  physical  forcing 
other  than  buoyancy,  pressure,  and  Coriolis.  The  following  discussion  thus  conveys  only  a  small  portion  of 
the  EU LAG’s  capabilities,  but  the  summarized  methodology  is  representative  of  the  entire  model.  Further, 
given  the  physical  scope  of  this  paper  we  use  a  concise  symbolic  operator-form  description  of  the  governing 
equations;  for  the  complete  tensorial  expositions  refer  to  Smolarkiewicz  and  Prusa  (2005)  and  references 
therein* 

3.2,  Analytic  formulation 

To  address  a  broad  class  of  flows  in  a  variety  of  domains  -  with,  optionally,  Diriehlet,  Neumann,  or  peri¬ 
odic  boundaries  in  each  direction  -  we  formulate  (and  solve)  the  governing  equations  in  a  transformed  domain 
with  time-dependent  curvilinear  coordinates 

(7,x)  =  (/,^(/,x)).  (I) 

The  key  assumptions  are  that  the  coordinates  (/,x)  of  the  physical  domain  are  orthogonal  and  stationary  -  in 
particular,  Cartesian  in  this  paper  -  and  that  the  transformed  horizontal  coordinates  [x,y)  are  independent  of 
the  vertical  coordinate  z.  Given  the  transformation  in  (1),  the  adiabtic  incompressible  Boussinesq  equations 
for  a  salty  water  can  be  compactly  written  as  follows: 


3  The  acronym  NFT  labels  a  class  of  seco nd-order-accu rate  two-lime- level  algorithms  buih  on  nonlinear  advcclion  techniques  that 
suppress/ red uce/coTitrol  numerical  oscillations  characteristic  of  higher-order  linear  schemes;  it  is  meant  to  distinguish  these  algorithms 
from  classical  ccntercd-in-iime-and-space  linear  methods. 
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V  ■  (pV)  =  0, 

(2) 

(3) 

(4) 

(5) 

where,  because  of  the  coordinate  transformation,  the  physical  and  geometrical  aspects  intertwine  each  other. 
Insofar  as  the  physics  is  concerned:  v  =  denotes  the  physical  velocity  vector;  p,  j  and  te  denote  density, 

specific  salinity,  and  a  density-normalized  pressure;  g  is  the  acceleration  of  gravity  (vector),  and  f  the  vector  of 
Coriolis  parameter.  Primes  denote  deviations  from  the  geostrophically-balanced  ambient  (alias,  environmen¬ 
tal)  state  pc;  and  pQ  is  the  Boussinesq  reference  density. 

The  geometryof  the  coordinates  in  (1)  enters  the  governing  equations  as  follows:  in  the  mass  continuity 
Eq.  (2),  p*  =  p0C7  with  G  denoting  the  Jacobian  of  the  transformation;  whereas  in  the  momentum  Eq<  (3), 
G  symbolizes  the  renormalized  Jacobi  matrix  of  the_transformation  coefficients  ~  ( dx/dx );  V'  =  0/8x-,  and 
the  total  derivative  is  given  by  d/d7  =  0/87  +  v*  *  V,  where  v*  =  dx/d7  =  *  is  the  contravariani  velocity . 
Appearing  in  the  continuity  (2)  and  density  (4)  equations  is  a  weighted  solenoidal  velocity 


v*  =  V 


8x 

dp 


(6) 


which  readily  follows  (Prusa  et  ah,  2001)  from  the  generic  (tensor  invariant)  form  of  the  mass  continuity 
equation 


G"’(^  +  V'(^9‘))  =0.  (7) 

Use  of  the  solenoidal  velocity  facilitates  the  solution  procedures  because  it  preserves  the  incompressible  char¬ 
acter  of  numerical  equations,  regardless  of  the  time-dependency  of  the  transformed  coordinates.  While  numer¬ 
ous  relationships  can  be  derived  that  express  any  velocity  (solenoidal,  contravariani,  covariant,  or  physical)  in 
terms  of  the  other,  in  either  transformed  or  physical  coordinate  system,  a  particularly  useful  transformation 

Vs  =  G 1  v,  (8) 

relates  the  solenoidal  and  physical  velocities  directly.  For  further  details  of  the  metric  and  transformation 
tensors  as  well  as  formulating  viscous  and  dissipative  terms  in  the  governing  equations,  the  interested  reader 
is  referred  to  Smolarkiewicz  and  Prusa  (2005)  and  the  references  therein. 

The  equations  of  motion  (2}-(5)  are  supplemented  with  a  linearized  constitutive  law  for  sea  water 

ff  =  « T  +  ps'  (9) 

with  constants  a  =  -2.5  x  10  4  and  f$  “  7,6  x  10  4,  representative  of  conditions  in  the  Mediterrean  region 
(Benoit,  1994;  Marshall  ct  aL,  1997).  The  linear  relation  (9)  is  used  to  diagnose  temperature/potcntial-temper- 
ature  perturbations  and  to  construct  the  ambient  density  profile  pe  from  the  potential  temperature  and  salinity 
profiles  and  % 

The  numerical  apparatus  employed  to  solve  the  posed  equations  is  summarized  in  Appendix.  Here  we  note 
that  there  are  four  important  benefits  of  formulating  governing  PDEs  as  described:  (i)  the  relative  simplicity  of 
designing  a  flow  solver  fully  implicit  with  respect  to  gravity  waves  (Appendix);  (ii)  conservation  of  density 
perturbations  with  accuracy  to  round-olT  error  -  cf.  Section  3a  in  Smolarkiewicz  et  al  (2001),  for  a  discussion 
-  tantamount  to  preventing  dilution  of  the  ambient  stratification  due  to  implicit  viscosity  of  nonoscillatory 
advection;  (iii)  improved  accuracy  of  impermeability  conditions  imposed  along  irregular  boundaries  (Smol¬ 
arkiewicz  et  aL,  2001);  and  (iv)  improved  energetics  of  the  model  (Wedi,  2006).  For  a  validation  against  an 
akin  laboratory  experiment  with  density  stratification  and  transient  nonlinear  gravity- wave  dynamics  see  Wedi 
and  Smolarkiewicz  (2006). 
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33.  The  model  setup;  details  of  transformation 

Following  Wed!  and  Smolarkiewicz  (2004),  the  general  dependence  of  z  on  {x,y,z, /)  in  (1)  collapses  into  a 
similarity  transformation 


{  =  :=  Ifo 


z  -  zt{x,y,  t) 
H(x,y,  t )-zs(x,y,t)  ’ 


(10) 


where  H  and  are  the  upper  and  lower  surface  elevations,  respectively,  //0  denotes  the  vertical  extent  of  the 
transformed  model  domain,  and  the  function  %  conveniently  admits  a  class  of  vertically  stretched  coordinates. 
The  transformation  in  (10)  is  a  generalization  of  the  classical  terrain-following  Gal-Chen  and  Somerville 
(1975)  transformation.  It  has  the  computational  advantage  of  separability  into  one-  and  two-dimensional 
fields.  In  particular,  the  Jacobian  of  the  transformation  is  given  as 


/ d#  dc\  1  /c \x  dy  dx  dy\  1 

Vz) 


with 

g  _  mV*  _  If{x,y,t)-z,(x.y,t) 
“  ~~  \02/  H0 


00 


(12) 


The  similarity  transformation  (10)  is  a  convenient  vehicle  for  bringing  a  varying  ocean  surface  from  other 
sources,  either  data  or  models.  For  illustration  see  Wedi  and  Smolarkiewicz  (2004),  where  the  authors  simu¬ 
lated  nonhydrostatic  flow  of  a  homogeneous  shallow'  layer  of  water  past  a  hill,  with  a  variable  upper  boundary 
H{xtytt)  in  (10)  predicted  by  integrating  the  shallow'-watcr  equations.  Furthermore,  the  dependence  of*  and  y 
on  the  horizontal  coordinates  of  the  physical  space,  admits  studies  of  fluid  flows  in  curved  channels  (Smol¬ 
arkiewicz  and  Prusa,  2005),  Notwithstanding  the  generality  of  (10),  throughout  this  paper  x  x, 
y  —  y  and  c  ~z;  thereby  employing  the  identity  transformation  in  the  horizontal  (viz.  Giy  =  1 ).  Furthermore, 
the  upper  boundary  is  stationary  and  flat  (viz,  H  =  Hq)<  and  there  is  no  vertical  stretching  of  the  lower-bound- 
ary-fitted  coordinate  z  (viz,  d#/de  =  1).  The  lower  boundary  is  also  stationary  but  inhomogeneous, 
rs  =  z$(x,y),  thereby  reducing  ( 10)  to  the  classical  case,  standard  in  many  atmospheric/oceanic  models.  In  spite 
of  the  resulting  mathematical  simplifications,  the  actual  EULAG  program  accommodates  (1)  and  (10)  in  their 
full  generality.  Consequently,  we  retain  the  consistent  notation  for  future  reference  and  ease  of  connection  to 
earlier  works. 


3,4.  The  model  setup;  ambient t  initial  and  boundary  conditions 

A  specific  perturbation  form  of  PDEs  solved  in  EULAG  depends  on  the  admitted  class  of  ambient  states. 
Consider  for  illustration  that  the  governing  system  (2)— (5)  derives  from  the  generic  Boussincsq  form  in  which 
the  pressure  and  temperature  perturbations  in  the  momentum  equation  are  taken  with  respect  to  the  static 
vertical  profiles  of  the  Boussincsq  expansion  -  simply  by  postulating  there  exists  an  inertia!  ambient  state 
determined  by  the  balance  of  pressure,  buoyancy,  and  Coriolis  forces 

0  =  —  G(V(jic  —  fin))  +  g£s_Z_fl£  -  f  x  vc,  (13) 

Pq 

together  with  implied  compatibility  conditions  (e.g.,  for  a  sheared  ambient  flow,  pe  must  change  accordingly  in 
the  horizontal;  cf.  Smolarkiewicz  et  ah,  2001,  for  a  discussion).  Subtracting  (13)  from  the  generic  form  results 
in  (3),  while  the  perturbational  formulation  of  the  density  Eq.  <4)  follows  readily  the  generic  form  dp/di  =  0, 
given  pe  =  pe(x). 
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This  feature  of  the  model  is  employed  to  incorporate  the  effects  of  tidal  forcing.  In  particular,  the  ambient 
flow  is  assumed  as  vc  =  [ue(r),0,0]f  where 


Here,  constant  VT  denotes  the  flow  magnitude  at  the  south  boundary  of  the  model  located  on  the  Ionian  side 
of  the  Messina  Strait,  with  numerical  value  (Table  1)  determined  after  Martin  (2000);  whereas  T  —  1 2,4  h  is  the 
semidiurnal  period  of  the  tidal  forcing.  Because  the  ambient  flow  is  allowed  to  vary  in  time,  the  governing 
momentum  equations  should  include  d\c/dl  on  the  rhs  -  in  light  of  the  preceding  discussion.  However,  since 
the  characteristic  time-scale  of  solitary- wave  evolution  is  much  shorter  than  F,  these  terms  are  neglected  in 
further  discussion.  The  perturbation  form  of  the  governing  equations  is  incapable  per  se  of  appreciating 
the  tidal  flow,  as  it  is  simply  a  result  of  subtracting  a  prescribed  subset  of  solutions  from  both  sides  of  the 
complete  Boussinesq  problem.  The  system  is  actually  connected  to  tidal  flow  via  inflow  boundary  conditions 
v  n  —  VT  imposed  at  the  south  inflow  boundary,  determining  there  the  Neumann  boundary  condition  for 
pressure  in  the  elliptic  boundary  value  problem  (Appendix).  In  order  to  account  for  the  sign-change  of  the 
tide  after  F/2,  at  the  time  of  the  model  initialization  the  constant  V  r  is  determined  at  the  north  lateral  model 
boundary,  such  as  to  satisfy  the  integrability  condition  fdQp*v*  ■  ndu  =  0  implied  by  (2).  As  the  tide  changes 
sign  and  the  south  and  north  boundaries  become  outflow  and  inflow,  respectively,  the  solution  accounts  for 
the  tide  via  v  n  =  FT  enforced  at  the  north  boundary.  To  minimize  numerical-boundary  artifacts,  dissipative 
absorbing  layers  attenuate  the  solution  gradually  towards  the  ambient  flow  in  the  vicinity  of  the  boundaries. 

To  ensure  a  well-posed  initial  condition  (cf.  Temam,  2006},  v{/  =  0)  is  assumed  a  superposition  of  the  ambi¬ 
ent  flow  and  potential  perturbation 

V(f  =  or=vc(f=o)-G(m  (is) 

and  it  is  determined  by  solving  (on  the  model  grid)  a  discrete  elliptic  problem,  for  the  potential  0,  implied  by 
the  mass  continuity  equation  (2)  and  the  imposed  boundary  conditions.  The  initial  fields  of  potential  temper¬ 
ature  and  salinity  are  set  to  ambient  values;  whereupon,  the  initial  density  perturbation  is  set  to  zero.  To  con¬ 
struct  the  initial  profiles,  the  ambient  temperature  and  salinity  fields  were  obtained  from  CTD  stations  shown 
in  Fig.  1.  For  modeling  convenience,  analytic  functions  were  fitted  to  the  measurements,  as  shown  Fig.  2. 

The  model  domain  of  Lx  x  Lv  x  H  =  80  x  10  x  0.4  km1  in  ,x,  y  and  r,  respectively,  is  covered  with 
800  x  100  x  100  grid  points.  The  vertical  extend  is  restricted  to  the  top  400  m  of  the  ocean.  The  ^-direction 
is  oriented  from  the  Strait  of  Messina  towards  the  Gulf  of  Gioia;  consult  Fig.  I  for  the  outline  of  the  horizon¬ 
tal  domain.  Numerical  simulations  cover  2-5  semidiurnal  tidal  periods  using  time  step  5 1  of  l  ,5  s  limited  by  the 
advective  CFL  condition.  The  z  =  zs  and  z—H  boundaries  arc  impermeable.  At  jc  —  0  the  lateral  boundary  is 
open,  whereas  at  x  =  Lx  the  bathymetry  zs(xt y)  continues  into  the  coastline.  To  prevent  a  singularity  in  the 
transformation  (10)  where  the  bathymetry  protrudes  to  the  ocean  surface  z  =  //,  the  actual  bathymetry  is  lim¬ 
ited  such  that  zs  =  mi n(z^H  -  A)  with  A  =  20  m;  cf.  Fig.  3.  Concomitantly,  in  the  sea  region  above  the 
maximum  bathymetry  zsmax  =  //  -  d,  flow  is  attenuated  rapidly  to  stagnation  using  the  Rayleigh  friction 
(on  the  rhs  of  the  momentum  equation)  -  in  the  spirit  of  gravity- wave  absorbers  common  in  atmospheric  mod¬ 
els  -  with  the  c-folding  time  scale  t  =  18  s,  At  the  y  —  Ly  model  boundary,  for  all  x  <  0.5 Lx  and  z  >  zs  the 
coastline  is  steep,  thus  justifying  the  approximation  with  impermeable  vertical  walls.  At  y  =  Lv  and  x  >  0..5Lr 
the  boundary  is  open  for  all  z  >  zs.  The  latter  results  in  a  complex  shape  of  the  opening  towards  Tyrrhenian 
Sea  -  an  intersection  of  the  (x,z)  plane  at  v  —  Lv  with  the  bathymetry.  This  motivates  the  customized  scheme 


Tabic  I 

Simulation  parameters:  the  cen  ter  of  the  thcrmocline  depth  Pjr  the  center  of  the  haloeline  depth  P$>  and  the  semidiurnal  tidal  magnitude 


Case 

fiN 

Mm] 

rT[ins-1] 

l 

40 

50 

0.5 

2 

40 

50 

0.25 

3 

60 

70 

0.5 

104 


A  Warn*  Varum  ei  at.  /  Ocean  Modelling  I  ft  (2007)  97-121 


a 


K4D. 

? 


HO. 


-40.  -24.  -B  8-  24.  40. 

x  (km) 


x  (km) 


*  Oan)  i  (km) 

Fig.  3.  Siraii  of  Messina  sill  with  Tyrrhenian  Sea  on  ihe  right  and  Ionian  Sea  on  the  left.  Isohalines  contours  from  37.995  ppl  to  38.638  ppi 
in  intervals  of  0.09  ppt  in  the  y  =  5  km  jr;  plane  at:  (a)  0.52  h;  (b)  1 .55  h;  (c)  3.62  h  and  (d)4.!3h. 


for  handling  the  tide  reversal,  discussed  earlier  in  this  subsection.  Aiy  —  0  the  boundary  is  impermeable  for  all 
r  >  rs,  and  the  narrowing  of  the  channel  at  the  sill  is  represented  by  the  bathymetry,  analogously  to  the  x  —  Lx 
boundary;  cf.  Fig.  7a. 

All  calculations  reported  in  this  paper  are  performed  using  the  Eulerian  flux-form  option  of  EULAG,  exe¬ 
cuted  in  the  ILES  mode,  in  which  the  truncation  terms  of  nonoscillatory  transport  schemes  play  the  sole  role 
of  the  subgrid-scale  turbulence  models;  see  Smolarkicwic2  and  Margolin  (2007)  for  a  review'. 


4.  Dynamics 

4  J.  Generation  of  interna!  bores,  fronts  and  solitary  waves 

At  the  Strait  of  Messina  silt  region,  the  salty  water  originates  from  Levantine  Intermediate  Water  (LIW), 
The  surface  water  can  consist  of  modified  allantic  w'ater  (MAW)  or  Tyrrhenian  Sea  surface  water  (TSW),  As 
the  semidiurnal  tide  moves  over  the  Strait  of  Messina  sill  from  the  Ionian  to  the  Tyrrhenian  Sea,  the  pycno- 
cline  and  halocline  region  is  moved  upwards  at  the  sill  and  depressed  on  the  right  side  of  the  sill.  Fig.  3a.  The 
depression  is  an  internal  bore  with  salinity  and  temperature  gradients  along  its  boundaries.  The  first  internal 
bore  to  be  formed  by  the  tidal  flow  to  the  right,  generates  a  right  and  a  left  propagating  internal  wave.  The 
salinity  field  for  such  an  internal  bore  is  shown  in  Fig*  3a  at  the  simulation  time  t  =  0.52  h  along  a  vertical  jcz 
plane  located  at  mid  distance  along  the  y  coordinate,  y  =  5  km.  The  parameters  of  the  simulation  are  listed  as 
case  1  in  Table  1,  where  PT  is  the  center  of  the  thermodine  depth,  Ps  is  the  center  of  the  halocline  depth,  and 
VT  is  the  prescribed  semidiurnal  tidal  magnitude  (cf  Section  3.4  for  discussion).  At  L55  h.  Fig.  3b,  the 
right  propagating  wave  steepens  on  the  leading  edge  through  nonlinear  effects.  Then*  later  in  time,  frequency 
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and  amplitude  dispersion  set  in  and  it  disintegrates  into  solitary  waves,  Figs.  3c  and  3d.  These  are  solitary 
waves  of  depressions  that  push  the  isohalines,  isotherms  and  isopycnals  down.  When  the  semidiurnal  tidal 
flow  reverses  the  direction,  from  positive  to  negative  along  x,  the  left  propagating  internal  bore  of  depression 
that  was  formed  by  the  positive  tidal  flow.  Fig.  3b,  undergoes  a  hydraulic  jump  over  the  sill.  Figs.  3c  and  3d. 

Inside  the  internal  bore.  Figs.  3c  and  3d,  there  is  a  tendency  for  a  more  homogeneous  temperature  and 
salinity  distribution.  At  the  bottom,  left  and  right  sides  of  it  there  are  large  salinity  and  temperature  gradients 
-  the  temperature  field  has  a  similar  distribution  to  the  salinity  field.  As  the  bore  moves  over  the  sill,  the  bot¬ 
tom  part  generates  a  front  between  the  surface  and  the  top  of  the  sill.  During  the  dynamics  of  the  motion  that 
takes  place  over  the  semidiurnal  tidal  cycle,  a  large  salinity  and  temperature  gradient  can  be  located  anywhere 
between  the  sill  and  the  surface  of  the  ocean.  The  water  masses  that  exist  in  the  lower  and  upper  parts  of  the 
sill  region,  are  separated  by  this  frontal  temperature  and  salinity  gradient.  At  t  —  6.2  h  the  salinity  distribu¬ 
tions  has  evolved  to  the  configuration  shown  in  Fig.  4a. 

There  is  a  salinity  gradient  or  front  located  past  the  mid  depth  over  the  sill.  The  two  water  masses,  upper 
and  lower,  are  on  each  side  of  the  front.  The  internal  bore  that  jumped  over  the  sill.  Fig.  3d,  has  now  moved 
close  to  the  x  —  -24  km  location  and  has  became  steep  on  the  leading  side  through  nonlinear  effects.  Proceed¬ 
ing  to  /  —  7,75  h  the  flow  over  the  sill  is  moving  to  the  left,  towards  the  Ionian  Sea,  As  the  tidal  flow  continues 
to  the  left,  the  salinity  front  moves  further  up  on  the  sill  towards  the  surface.  The  gradient,  associated  with  the 
front,  is  now  located  near  the  surface,  Fig.  4b.  To  the  right  of  the  sill,  the  solitary  wave  train  that  was  previ¬ 
ously  near  x  =  24  km,  Fig.  4a,  is  heading  towards  the  Gulf  of  Gioia,  Fig.  4b. 

The  location  of  the  salinity  front  changes  as  one  proceeds  westward  from  the  mid  plane  {at  >'  =  5km) 
towards  increasing  y  values.  In  a  vertical  plane  close  to  the  Sicilian  shore,  y  —  8  km,  at  the  time  of  7,75  h 
the  front  is  anchored  on  the  left  side  of  the  sill  and  extends  to  the  surface  with  a  upwards  slope.  Fig.  5a.  When 
the  tidal  flow  reverses  to  be  along  the  positive  ^-direction,  at  /  =  1 1.37  h,  the  front  becomes  anchored  on  the 
right  side  of  the  sill  and  extends  to  the  surface  by  sloping  upwards  to  the  left.  Fig.  5b.  The  low  salinity  water 
tends  to  be  on  the  right  side  of  the  front  and  the  high  salinity  water  on  the  left.  In  the  previous  case.  Fig.  5a, 
the  reverse  situation  existed.  The  solitary  wave  train  that  was  located  at  x  =  24  km  in  Fig,  5a  at  t  —  7.75  h,  has 
reached  the  Gulf  of  Gioia  shore  at  t  =  1 1.37  h,  shoaled,  and  resulted  in  a  depression  along  the  shelfbreak. 
Fig.  5b, 

Hopkins  et  al,  (1984)  and  Del  Rieco  (1981)  have  investigated  the  time  evolution  of  the  interface  between  the 
two  water  masses  over  the  Messina  sill  with  two-layer  models.  The  interface  of  the  two-layer  models  was 
found  to  exhibit  large  vertical  oscillations  between  the  surface  and  bottom  of  the  sill  region.  The  interface 
can  be  at  the  surface,  bottom  or  anywhere  in-between.  Our  simulation  with  the  nonhydrostatic  model  exhibit 
a  similar  phenomena  of  water  mass  separation  by  a  frontal  gradient  of  temperature  and  salinity  caused  by  the 
semidiurnal  tidal  motion  over  the  sill. 

The  solitary  waves  trains  propagate  away  from  the  sill  The  structure  of  the  solitary  waves  changes  along 
and  across  the  direction  of  propagation.  Consider  a  vertical  section  located  towards  the  outer  boundary  in  the 


Fig.  4  Isohalines  contours  from  37.975  ppl  to  38.636  ppt  in  intervals  of  0.094  ppi  in  the  y  =  5  km  xs  plane  at:  (a)  6.2  h  and  (b)  7.75  h. 
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Fig.  5.  Iso  ha  lines  contours  from  37.993  ppt  to  38.642  ppt  in  intervals  of  0,092  ppt  in  a  at  plane  at  r  =  8  km  at:  (a)  t  =  7.75  h  and  (b| 
i =  1 1J7  h. 

y-direction,  Fig.  5a,  The  beginning  of  the  solitary  wave  train  at  /  =  7.75  h  is  near  ,v  =  24  km.  Along  the  mid 
section  of  the  v-dtrectiom  Fig.  4b,  the  solitary  wave  train  is  located  past  x  =  24.0  km,  which  is  further  ahead  in 
the  A-direction.  As  one  moves  across  the  direction  of  propagation  in  the  ^-direction,  there  is  a  curvature  effect 
in  the  location  of  the  solitary  wave  trains.  The  curvature  of  the  solitary  wave  train  can  be  seen  in  a  horizontal 
plane  located  100  m  below  the  surface.  Fig.  6.  At  around  65  km  down  range  there  is  a  warmer  temperature 
distribution  generated  by  the  depressions  associated  with  the  solitary  waves.  These  depressions  bring  the  war¬ 
mer  surface  temperature  down.  Note  the  curvature  of  the  warmer  temperature  arcs.  The  arcs  curve  backwards, 
as  one  proceeds  from  land  (bottom)  towards  the  open  sea  (top), 

42.  Flow  structure 

At  /  —  7.75  h  consider  the  semidiurnal  tidal  circulation  that  exists  in  the  numerical  domain.  The  forcing  of 
the  semidiurnal  harotropic  flow,  Eq,  (14)  in  Section  3.4,  is  in  the  reverse  direction  at  this  lime.  In  the  sill  region 
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Fig.  6.  Horizontal  temperature  distribution  at  7.75  h  in  a  plane  located  100  m  below  ihe  surface.  Down  range  corresponds  to  the  x- 
direction  and  cross-range  to  the  v-dircetion. 
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Fig.  7r  I  lorizontal  components  or  velocity  at  t  —  7r75  h:  {a)  w  component  in  the  yz  plane  at  30  km  down  range  of  Fig.  6  with  contours  from 
-2.75  to  “0.33  ms'1  in  intervals  of  0.04  ms_I;(b)it  component  at  the  60  km  down  range  with  contours  from  -0.363  to  +0.53  m  s  1  in 
intervals  of  0.128  ms-1  and  (c)  y  component  in  the  same  plane  as  (b)  with  contours  from  -0.203  to  -0.013  ms“!  in  intervals  of 
0.047  ms*1.  Vertical  and  horizontal  axes  are  marked  in  m  and  km.  respectively. 


between  Italy  and  Sicily,  Fig.  7a  at  30  km  down  range  (over  the  sill,  cf.  Fig.  6)  there  is  a  large  section  of 
negative  flow  (towards  the  Ionian  Sea)  of  around  3ms  1 .  The  yellow  upper  corners  in  the  figure  arc  modeling 
artifacts  that  represent  continuation  of  the  bathymetry  by  means  of  stagnant  water;  cf.  discussion  at  the  end  of 
Section  3  A  Further  down  the  x-axis  at  60  km,  between  the  sill  and  the  Gulf  of  Gioia,  the  flow  structure  has  a 
large  section  directed  towards  the  sill,  Fig.  7b,  extending  from  the  right  to  left  boundaries  and  encompassing 
the  bottom.  The  flow  towards  the  sill  on  the  right  boundary  reflects  the  inflow  from  the  side  boundary  con¬ 
dition.  On  the  left  boundary  there  is  flow  towards  the  sill,  also,  in  the  lower  layer.  This  indicates  a  local  cir¬ 
culation,  In  the  yz  plane  at  40  km  down  range  (not  shown),  the  flow  away  from  the  sill  (in  positive  x  direction) 
disappears  and  all  the  flow  is  towards  the  sill. 

The  flow  is  superposed  on  a  sloping  bottom  that  deepens  as  one  proceeds  away  from  the  land  into  the  sea. 
At  the  60  km  down  range.  Fig.  7c,  the  spanwise  velocity  indicates  inflow  along  the  ocean  y  =  Ly  boundary 
where  the  barotropic  forcing  flow  is  prescribed  {Section  3,4).  The  spanwise  flowr  controls  the  negative  flow 
away  from  the  sill  in  the  x-direction,  When  the  semidiurnal  tidal  forcing  reverses,  the  direction  of  flow  on 
the  sill  and  other  locations  reverses  also. 


43.  Shoaling  in  the  Gulf  of  Gioia 

The  solitary  wave  trains  that  are  generated  in  the  Strait  of  Messina  region  propagate  into  the  Gulf  of  Gioia. 
Fig.  4  shows  a  solitary  wave  train,  in  the  Gulf  of  Gioia  area,  approaching  the  shore  over  a  rising  topography. 
The  fiat  region  at  the  right  represents  continuation  of  the  bathymetry  by  means  of  stagnant  water;  Section  3.4. 
At  a  time  of  8,78  h  the  train  has  moved  closer  to  shore.  Fig.  8a.  The  initial  soliton  in  the  train  has  been  slowed 
down  by  the  topographic  obstruction  and  the  subsequent  soli  tons  are  catching  up  with  it.  Eventually  the  sol- 
itons  in  the  train  are  pushed  together  and  one  resultant  depression  is  formed  along  the  topography,  Fig.  8b, 
This  depression  is  slanted  along  the  topographic  rise.  There  is  a  pushing  down  of  less  dense  water  to  depths  of 
around  150  m  from  the  surface  over  a  horizontal  extend  of  about  1^2  km.  This  opens  the  possibility  for  water 
mass  modification  through  mixing.  At  a  later  time  of  13.95  h  the  depression  along  the  topography  disappears 
and  the  isopycnals  move  up,  towards  their  original  equilibrium  position,  Fig.  8c,  There  remains  a  signature  on 
the  isopycnals  to  the  left  of  the  shore  line  suggesting  reflectance  effects.  Reflectance  effects  have  been  observed 
by  Guardian!  et  ai.  (1988). 

In  a  horizontal  plane  at  a  depth  of  100  m,  the  temperature  field  indicates  the  disturbances  induced  by  the 
solitary  waves  and  internal  bores,  Fig.  9.  The  plate  (a)  corresponds  to  Fig.  8a  that  shows  the  solitary  wave 
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Fig.  8.  Shoaling  of  a  solitary  wave  train,  Isohatinc  contours  from  37.992  ppi  to  3K.63H  ppt  in  intervals  of  0.092  ppt  in  the  v  =  5  km  j.vr) 
plane  at;  (a)  i  =  8,78  h;  (b)  i  =  10.33  h  and  (c)  t  =  1 3*95  h„ 


train  about  to  shoal.  Fig.  9a  exhibits  warmer  temperatures  caused  by  the  solitary  wave  train  depressions  and  a 
curvature  of  the  solitary  waves  away  from  land,  bottom  of  graph.  At  a  time  of  1 0 33  h,  Fig.  9b,  the  depressions 
of  the  solitary  waves  line  up  along  the  shelfbreak  slope  and  result  in  a  warmer  temperature  there.  At 
/  —  13.95  h,  Fig,  9c,  the  depressions  along  the  shelfbreak  slope  have  disappeared  and  the  temperature  distri¬ 
bution  has  returned  towards  its  normal  values. 

We  have  introduced  three  moorings  into  the  domain  of  the  numerical  experiment.  The  moorings  were 
placed  at  63,  68  and  69,5  km  down  range,  in  the  middle  (xz)  plane  at  5  km  cross-range.  The  mooring  at 
68  km  is  centered  in  the  leading  depression  of  the  shoaling  solitary  wave  train  shown  in  Fig,  8a  at 
f  =  8.7'8h.  The  resultant  temperature,  salinity,  and  density  distributions,  at  the  moorings,  are  shown  in 
Fig.  10  at  the  times  of  8.78  h,  10,33  h,  and  13.95  h.  The  density  distribution  for  t  =  8.78  h  indicates  a  down¬ 
ward  displacement  of  the  pycnocline  for  the  mooring  at  68  km  dowm  range,  green  curve,  relative  to  the  moor¬ 
ing  at  63  km,  blue  curve.  For  i  =  10.33  h,  the  lime  at  which  a  large  depression  along  the  shelfbreak  slope  is 
evidenced  in  Fig.  Sb,  there  are  lighter  density  regions  at  depth  for  the  mooring  location  of  68  and  69.5  km 
down  range,  green  and  red  curves.  This  suggests  mixing  below'  the  surface.  At  /  -  13.95  h  conditions  relax 
towards  normal. 

Previously,  Guardian!  et  a L  (1988)  have  noted  a  patch  of  mixed  water  at  mid  depth  in  the  Gulf  of  Gioia, 
during  the  JANE  1984  cruise.  The  patch  was  centered  from  40  m  to  100  m  in  depth.  It  was  suggested  that  tur¬ 
bulent  mixing  due  to  breaking  of  internal  waves  could  be  the  cause  of  a  patch  of  mixed  water.  During  the 
shoaling  of  internal  wave  trains,  we  do  not  observe  internal  wave  breaking  in  the  model  predictions.  The 
shoaling  consist  of  solitary  weaves  of  the  train  piling  up  along  an  extended  shelf  slope  causing  a  large  depres¬ 
sion  that  eventually  return  to  equilibrium.  What  we  do  observe,  however  is  a  movement  down  of  the  surface 
water  mass  and  a  subsequent  movement  back  up  to  the  original  position.  This  up  and  down  movement  causes 
water  mass  modifications  through  mixing,  at  depths  of  about  30  m  to  150  m. 
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Fig.  9,  Horizontal  temperature  distribution  in  a  plane  located  1 00  m  below  the  surface  at:  (a)  t—  8.78  h;  (b)  /=  10,33  h  and  (c) 
t  —  13.95  h.  Down  range  corresponds  to  the  .v-di  reel  ion  and  cross-range  to  the  v-direction. 
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Fig.  10.  Model  predicted  tempera  lure,  salinity,  and  density  profiles  Tor  moorings  located  in  the  central  Ur)  plane  at  down  range  63  km 
{blue  curve),  68  km  (green  curve)  and  69.5  km  [red  curve)  for  times  r  =  8.78  h,  /  -  10.33  h.  and  t  =  13.95  h. 
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5*  Comparisons  with  data  and  parameter  sensitivities 

5.  /.  Wavelet  analysis  of  soli  ton  packets 

Traditionally  researchers  have  relied  on  standard  Fourier  techniques  to  investigate  the  spectral  composition 
of  soliton  packets.  Nonetheless,  because  Fourier  analysis  gives  signal  amplitudes  averaged  over  entire  tempo¬ 
ral  or  spatial  domain,  it  is  mdiscrimmating  when  a  signal  evolution  is  concerned.  Unlike  global  Fourier 
decomposition,  wavelet  analysis  expands  one-dimensional  time/space  series  into  the  twrQ-dimensional  param¬ 
eter  space  (a, 6)  to  provide  a  local  measure  of  relative  amplitude  of  activity  at  scale  a  and  temporal  or  spatial 
location  h  {Meyers  et  al.,  1993;  Torrence  and  Compo,  1998).  In  order  to  obtain  the  characteristic  scales  and 
phase  speeds  exhibited  by  simulated  soliton  packets,  we  applied  the  wavelet  analysis  to  isopycnal  displace¬ 
ments  in  the  pycnocline.  We  use  the  Morlct  wavelet,  proven  successful  in  analyzing  dispersion  of  Yanai  waves 
in  an  ocean  model  (Meyers  et  al.,  1993). 

The  analyzed  signal,  Fig.  1 1,  represents  a  simulated  soliton  train  in  range  (x-direction)  consisting  of  1024 
grid  points  linearly  interpolated  from  the  model  output  to  obtain  a  resolution  of  10  m.  The  scale  parameter  a  is 
varied  to  obtain  wavelengths  extending  from  103  m  to  2.81  km.  Because  the  increase  in  sampling  resolution  is 
an  artifact  of  the  interpolation,  the  scales  considered  in  practice  range  from  a  few  100  m  to  just  over  2  km. 

Consider  the  mid-thermoclinc  displacement  at  /  =  6,2  h  in  Fig.  11a.  This  is  identically  the  result  shown  in 
Fig,  4a,  however  the  range  has  been  shifted  in  order  to  measure  the  solitary  wave  train  from  the  sill  and  facil¬ 
itate  comparison  to  data.  The  contour  indicates  a  train  with  several  depressions.  The  wavelet  power  spectra  of 
this  train  is  shown  in  Fig.  1  lb.  The  figure  exhibits  the  power  as  a  function  of  range  and  wavelength.  The  abso¬ 
lute  maximum  intensity  occurs  at  a  wavelength  near  1.33  km,  the  range  location  corresponding  to  the  back 
side  of  the  first  depression  at  31  km.2  A  tracking  of  the  maximum  power  intensity  as  a  function  of  wavelength 
is  shown  by  the  black  curve  in  Fig.  11c,  also  referred  to  as  the  spine.  It  exhibits  wavelengths  ranging  from 
about  0.75  km  to  2  km,  where  intensities  below  the  background  intensity  of  749.88  are  not  shown.  The  var¬ 
iation  of  the  wavelengths  along  the  spine  over  the  contoured  area  show  initial  decrease  of  location  in  range 
as  the  small  wavelengths  increase  to  larger  ones,  Fig.  1  lc.  This  section  of  the  spine  contains  large  wavelength 
gradients.  The  variation  reaches  a  minimum  value  in  range,  along  the  spine.  The  location  of  the  minimum 
value  in  range  is  in  the  vicinity  of  large  wavelength  gradients,  Fig.  1  lc.  After  the  minimum,  the  range  increases 
along  the  spine  as  the  wavelengths  increase. 

The  behavior  of  the  spines  at  t  =  4J3,  6.2,  and  7,75  h  are  shown  in  Fig.  12b  together  with  the  model  pre¬ 
dicted  soliton  train  structures  in  the  pycnocline,  Fig,  12a.  The  spines  are  shown  as  a  function  of  wavelengths 
versus  range;  the  range  spans  20-45  km  in  order  to  show  the  evolving  spine  pattern.  Each  of  the  spines  has  a 
minimum  range  value.  The  wavelengths  derived  from  the  maximum  wavelet  power  spectra  vary  as  1.16,  1.33, 
and  1.4  km,  Table  2.  The  wavelength  represents  distances  between  troughs  of  the  solitons,  As  time  progresses 
the  distance  between  the  first  two  troughs  of  the  soliton  train  tends  to  increase.  Fig.  12a.  The  amplitude  of  the 
first  soliton  decreases  from  52.3  m  to  48.5  m  from  the  first  to  the  second  location  of  the  train.  This  follows  the 
expected  trend  of  decreasing  amplitude  in  range.  Then  it  increases  somewhat,  to  50.1  m,  at  the  third  location 
of  the  solitary  wave  train.  This  increase  can  be  a  fluctuation  related  to  a  decreasing  ocean  depth. 

Along  each  spine  the  wavelengths  increase  from  small  to  larger  ones.  High  wavelength  gradients  occur 
along  the  spine  slope  as  the  wavelength  increases  up  to  minimum  range  value  of  the  spine.  Fig.  1  lc.  The  loca¬ 
tion  of  this  section  of  the  spine  tends  to  extend  to  larger  wavelengths  as  time  progresses.  Fig.  12b.  Beyond  the 
minimum  range  point  along  the  spines,  there  are  two  dominant  spine  slopes  with  a  transition  region  between 
them,  Fig.  12b.  The  structure  of  the  power  distribution  contours,  Fig,  lie,  also  undergo  changes  as  time 
progresses.  The  elliptical-like  shape  of  power  contour  tends  (not  shown)  to  become  more  compressed  along 
the  wavelength  distribution  axis  and  more  elongated  in  range  as  time  increases.  The  major  axis  of  the  ellipse 
undergoes  rotation  and  deformation.  The  deformation  reflects  the  evolution  of  the  soliton  trains  in  time. 


2  Note  that  the  range  shift  results  in  a  roughly  1 1  km  difference  between  Figs.  1 1  and  4a,  where  the  first  depression  is  located  at  about 
20  km. 
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Fig.  11.  (a)  Model  predicted  mid-thcrmoclinc  displacement  at  lime  of  6.2  h:  (b)  displacement  intensity  [nr]  as  a  function  of  wavelength 
and  range  and  (c)  spine  (black  line)  i.e.,  maximum  intensity  for  each  wavelength  together  with  the  intensity  contours,  which  range  from 
749.88  to  1 1 ,248  and  span  I  L76  dB. 
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Fig,  12,  Soliton  evolution:  (a)  soliton  trains  as  functions  of  depth  and  range  at  f  ~  4. 1 3  h  (blue),  t  =  62  h  (green)  and  t  =  7.75  h  (red) 
using  the  mid-thcrmocline  displacement  and  (b)  the  corresponding  spines  as  a  function  of  wavelength  and  range. 


Table  2 

Model  predictions  for  case  1  in  Table  l;  kp  is  the  wavelength  corresponding  to  the  maximum  wavelet  power  spectra,  A\  is  the  amplitude  of 
the  1st  soliton  and  r  is  the  range 

Wavelength  2P  [km]  1.16  1,33  1.4 

Amplitude  A 1  [m]  52.3  48.5  50.1 

Range  r  [km]  24.78  30,96  36,08 


5.2.  Comparisons  with  data  characteristics 

The  evolution  of  the  measured  soliton  wave  trains  over  time  and  distance  are  compared  with  the  evolution 
of  the  model  predicted  solitary  wave  trains.  The  CTD  chain  data  described  in  Section  2  is  used  for  comparing 
the  model  results.  The  range  of  the  EULAG  model  predicted  solitary  wavetrains,  Fig.  12a,  are  different  from 
the  ranges  of  the  measured  solitary  wavelrains,  Fig.  1 3a.  The  differences  in  ranges,  or  distances  from  the  sill,  of 
the  model  results  and  the  measurements  are  due  to  the  tidal  phase  used  in  initialization  of  the  model  results. 
The  focus  of  the  comparison  is  not  on  the  range  location  but  on  changes  over  distance  and  time. 

For  the  three  CTD  chain  passes  on  the  Tyrrhenian  side,  the  inid-pycnocline  displacement  was  extracted  for 
the  wavelet  analysis.  Table  3  lists  the  wavelengths,  obtained  through  wavelet  analysis,  at  maximum  power  for 
tows  5,  6,  and  8  at  ranges  of  12.68,  15.65,  and  20.67  km.  The  maximum  power  wavelength  decreases  in  range 
and  then  increases.  This  is  different  from  the  model  behavior.  Tabic  2,  where  the  wavelengths  increase  with 


A  Warn-  Vamus  et  al  /  Ocean  Modelling  IS  (2007)  97  12! 


m 


Fig,  13.  (a}  Measured  soliton  trains  at  /  =  61 6,83  h  [blue,  tow  5),  /  =  617,95  h  (green,  tow  6),  and  t  =  619.47  h  (red.  low  8),  and  (hi 
corresponding  spines  as  a  function  of  wa%relengths  and  range. 


range.  We  attribute  the  behavior  of  wavelengths  extracted  from  data  to  the  presence  of  features  in  the  ocean 
such  as  eddies  that  can  compress  and  elongate  the  solitons. 

The  amplitude  of  the  first  solilon  decreases  from  the  first  soliton  train  to  the  second.  Table  3.  Front  the 
second  to  the  third  train,  the  situation  becomes  more  complicated.  The  first  soliton's  location  and  amplitude 
are  influenced  by  a  "feature”,  Fig,  13a,  The  location  of  the  first  soliton  is  advanced  relative  to  the  second  and 
its  amplitude  increases  relative  to  the  second  train.  The  magnitudes  of  the  measured  first  soliton  amplitudes. 
Table  3,  are  lower  than  the  model  predicted  amplitudes.  Table  2.  For  the  first  and  third  soliton  trains  the  dif¬ 
ference  is  around  6  m,  with  the  model  results  being  larger.  For  the  second  train  the  model  results  arc  about 
15  m  larger.  This  suggests  some  feature  distortion,  in  the  data  for  the  amplitude  of  the  first  soliton  in  the  sec¬ 
ond  train. 

The  behavior  of  the  spines  as  a  function  of  range  for  each  of  the  tows  (5,  6  and  8)  is  shown  in  Fig.  13b 
together  with  the  mid-pycnocline  displacement,  Fig,  13a.  The  spine  for  tow  number  6,  is  reminiscent  of  the 
spine  for  tow  number  5  with  the  modification  of  the  initial  minimum  in  range.  Wavelengths  below  0,5  km 
and  above  2.5  km  are  in  the  power  intensity  background.  For  tow  number  8,  the  spine  has  a  discontinuous 
decrease  in  range  for  wavelengths  ranging  from  about  0.7  km  to  L2  km.  Fig.  13b.  This  discontinuous  jump 
is  associated  with  the  solitary  wave  signal  distortion  that  is  visible  in  Fig.  13a,  where  the  distance  between 
the  first  and  second  depression  is  increased  relative  to  tows  number  5  and  6.  This  is  caused  by  a  feature. 

The  wavelengths  tracked  by  the  spines  arc  measured  by  the  tows  and  predicted  by  the  model  at  certain 
times.  A  curve  of  wavelength  as  function  of  range  and  time  can  be  obtained  from  the  spine  distributions.  From 
the  slope  of  these  curves,  the  phase  speeds  are  calculated  as  function  of  wavelength,  solid  blue  curve  in  Fig.  14. 
The  discontinuity  in  phase  speed  between  the  wavelengths  of  about  0,65  km  and  1.1  km  is  evident  and  involves 


A ,  Warn-Varnas  el  all  Ocean  Modelling  18  (2007)  97-121 


1 15 


Table  3 

Results  from  data:  /.r  is  the  wavelength  corresponding  to  the  maximum  wavelet  power  spectra,  A\  is  the  amplitude  of  the  1st  soliton  and  r 
is  the  range 

Wa vele  ngt  h  Ap  [km  ]  1 .09  0. 77  1  ,  68 

Amplitude  A\  [m]  46,58  33.3  44.57 

Range  r  [km]  12.68  15.65  20.67 


Fig.  14.  Calculated  phase  speeds  as  a  Function  of  wavelength:  green  curve  is  derived  From  model  predictions  and  blue  curve  From 
measurements. 


a  decrease  in  phase  speed  of  about  0.25  m  s  '.  This  discontinuity  can  be  caused  by  a  feature  that  results  in  the 
advection,  or  stretching,  of  the  first  depression  of  the  train  away  from  the  second  depression.  Fig.  13a.  The 
other  depressions  in  the  soliton  train  can  also  be  affected  by  the  feature.  Such  a  feature  can  be  a  local  eddy 
with  a  velocity  component  along  the  solitary  wave  train.  A  correction  for  the  movement  due  to  advection 
can  be  introduced  by  moving  the  discontinuous  portion  of  the  phase  speed  curve  up,  dotted  blue  line  in 
Fig,  14. 

The  model  predicted  phase  speed,  green  curve  in  Fig,  14,  has  similar  trends  to  the  corrected  measured  phase 
speed  curve.  From  wavelengths  of  0.7  km  to  1.6  km,  the  phase  speeds  increase  up  to  a  maximum,  decrease, 
and  then  increase  again.  In  the  comparative  range,  the  model  phase  speeds  range  from  0,93  m  s  3  to 
0.86  m  s  1  and  phase  speeds  derived  from  data  range  from  0,88  m  s” 1  to  0.77  ms  '.  For  the  comparative  por¬ 
tion  of  the  curve,  the  model  results  have  somewhat  higher  phase  speeds  and  a  slightly  different  dispersion  ver¬ 
sus  wavelength  variation.  The  overall  phase  speed  trends  are  similar  between  model  and  data.  For  the  non¬ 
comparative  portions  of  the  curves  -  below  wavelengths  of  0.7  km  and  above  1.6  km  -  the  phase  speed  trends 
of  model  predictions  and  data  are  different.  These  are  the  portions  where  the  power  intensity  distribution  in 
wavelength  and  range  merges  with  the  background;  of.  Fig.  1 1  b. 

Sapia  and  Salusti  (1987)  observed  internal  solitary  waves  in  the  Strait  of  Messina  and  derived  theoretical 
estimates  of  linear  and  nonlinear  phase  speeds  from  hydrography  measurements.  On  the  north  side  of  the 
strait  the  linear  phase  speed  estimate  was  0.6  ms  1  and  the  nonlinear  0,72  ms"1.  These  calculations  are 
dependent  on  the  measured  density  and  the  application  of  a  two-layer  approximation.  Measurement  with  a 
towed  temperature  sensor  yielded  nonlinear  phase  speed  estimates  of  0.8  ms-1.  Phase  speeds  of  0.8  ms  1 
are  in  the  range  of  the  phase  speeds  derived  from  the  present  data  sets. 

5.3.  Parameter  variations 

We  have  varied  the  initial  conditions  and  the  magnitude  of  the  semidiurnal  barotropic  tidal  forcing.  Table  1 
lists  the  parameters  variations  that  we  considered.  Left  and  right  plates  in  Fig.  2  illustrate  the  initial  temper¬ 
ature  and  salinity  profiles,  respectively,  that  were  used  for  model  initializations. 

A  comparison  of  derived  phase  speed  as  a  function  of  wavelength  for  different  parameters  is  shown  in 
Fig.  15.  Cases  1  and  the  data  are  the  same  as  in  Fig.  14.  Case  2  represents  a  lower  barotropic  tidal  forcing. 
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Fig.  15.  Phase  speed  as  a  function  of  wavelength  for  different  parameters.  Table  1  lists  (he  model  parameters  for  the  cases.  Case  I  and  data 
phase  speeds  are  as  m  Fig.  14r 

The  predicted  phase  speeds  versus  wavelength  are  smaller  than  those  for  case  1.  A  lower  tidal  forcing,  results 
in  less  isopycnal  displacement,  retardation  in  solitary  wave  formation,  and  smaller  tidal  background  velocity 
(Lamb,  1994;  Wam-Varnas  et  ai.,  2003).  As  a  result  the  phase  velocities  are  smaller.  Case  3  represents  deeper 
thermoclme  and  haloclines.  The  phase  speed  has  an  increasing  trend  and  for  wavelengths  exceeding  l  km  is 
larger  than  those  of  case  L  For  wavelengths  from  about  0.65  km  to  1  km,  the  phase  speed  is  smaller  than 
in  case  1  and  closer  to  the  data.  Fig.  15.  The  phase  speed  variation  versus  wavelength  exhibits  less  variation 
than  case  1  or  the  data.  Case  1  and  the  data  show  an  increase,  a  decrease,  and  an  increase  again  to  about  wave¬ 
lengths  of  1.6  km.  Case  I  exhibits  more  of  the  data  trends  than  case  2  and  3. 

The  wavelengths  derived  from  wavelet  analysis  are  summarized  in  Table  4  for  the  three  cases.  Case  1  indi¬ 
cates  wavelengths  ranging  from  1.16  km  to  1 .4  km.  Case  3  wavelengths  range  from  l  .25  km  to  2.18  km.  Case  3 
wavelengths  are  correspondingly  larger  for  each  of  the  three  solitary  wave  trains.  The  deeper  ihermocline  and 
halocline  locations  of  case  3  lead  to  larger  wavelengths  in  relation  to  case  L  The  first  soliton  amplitudes  for 
case  3,  Table  4,  are  larger  for  the  second  and  third  solitary  wave  trains  than  the  corresponding  amplitudes  for 
case  1.  For  the  first  soliton  train  the  amplitude  is  smaller  for  case  3  than  in  case  L  The  first  soliton  could  still 
be  in  the  development  stage. 

In  Case  2,  the  lower  tidal  forcing  at  the  boundary  results  in  decreasing  wavelengths  from  the  first  to  the 
third  solitary  wave  trains.  Table  4.  A  lower  tidal  forcing  results  in  smaller  internal  bores  that  disintegrate  into 
wider  width  solitons,  The  corresponding  amplitudes  of  the  first  soliton  increase  for  the  first  train  to  the  second 
and  then  decrease  in  the  third  train.  The  magnitude  of  the  amplitudes  are  smaller  than  those  for  case  1,  ease  2 
and  the  data.  A  smaller  tidal  velocity  leads  to  smaller  displacements  of  isopycnals  and  smaller  amplitudes  for 
the  solitons.  The  tidal  forcing  magnitude  of  case  2  does  not  lead  to  amplitudes  close  to  those  observed  in  the 
data. 
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Table  4 

Model  prediction  summary:  k  ihe  maximum  power  spectra  wavelength,  A\  is  the  amplitude  of  the  1st  soli  ton  and  r  is  the  range 


Case 

/.p  [km] 

Ax  [ra] 

r  [km] 

1 

1.16 

52.3 

24.78 

1.33 

48.5 

30.96 

1.4 

50.  t 

36.08 

2 

2.18 

17.49 

20.40 

2.05 

23.87 

25,71 

1.9 

21.32 

29.78 

3 

1.25 

439 

25.65 

1,67 

52.1 

33.1 

2.18 

51.27 

38.56 

6.  Acoustical  effects 

From  an  acoustic  point  of  view  the  interesting  pari  of  this  study  was  the  3D  aspects  of  the  modeled  ocean 
environment.  The  ocean  model  results  covered  an  area  of  10  km  by  80  km.  This  area  is  shown  by  the  dashed 
box  of  Fig.  1, 

Acoustic  model  simulations  were  used  to  illustrate  the  effects  of  a  3D  ocean  environment  as  compared  to 
the  acoustic  effects  of  a  2D  environment.  The  acoustic  simulations  were  made  using  the  finite  element  para¬ 
bolic  equation  (FEPE)  model  (Collins,  1988a;  Collins,  1988b;  Collins  and  Westwood,  1991).  The  FEPE  model 
gives  accurate,  wave- theoretic,  numerical  predictions  of  acoustic  pressure  at  ail  forward  propagating  angles 
about  the  direction  of  propagation.  Since  the  model  uses  a  marching  solution  technique  taken  at  intervals  that 
arc  a  fraction  of  the  acoustic  wavelength,  the  model  is  ideal  for  accurately  including  small  ocean  environmen¬ 
tal  changes,  and  their  effects  (refraction,  diffraction,  interference)  on  the  acoustic  field.  The  FEPE  model  has 
been  verified  as  benchmark  accurate  (Collins,  1990);  thus,  it  was  an  appropriate  choice  for  these  simulations. 

Three  tracks  were  chosen  that  would  illustrate  the  acoustic  effects  of  the  simulated  3D  ocean  environment 
in  Fig.  3d.  One  track  was  in  the  middle  of  the  ocean  area  (designated  as  “Center”).  One  track  was  parallel  to 
the  Center  track  and  displaced  4  km  to  the  west  (designated  as  “West”)  of  the  Center  track.  One  track  was 
parallel  to  the  center  and  displaced  4  km  to  the  east  (designated  as  “East”).  These  three  tracks  traversed  par¬ 
allel,  adjacent  ocean  environments  that  were  oriented  from  the  northeast  to  southwest  (refer  to  Fig.  1).  The 


Fig.  16,  Transmission  loss  vs.  range  plot  for  the  Center  (gray  line),  West  (dashed  line)  and  East  (solid  black)  tracks  for  a  frequency  of 
400  Hz,  source  depth  50  m  and  receiver  depth  of  90  m. 
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Range  (km) 

Fig.  17.  Transmission  loss  vs.  range  for  the  same  three  tracks  shown  in  Fig.  16,  but  with  ihe  bathymetry  the  same  for  all  three  tracks. 


starting  location  is  in  the  Gulf  of  Gioia  and  going  southwest.  A  source  depth  of  50  m  was  selected  for  each 
track  and  the  start  range  for  each  track  made  to  correspond  to  with  the  70  m  depth  contour  line. 

Fig.  16  shows  the  results  of  the  acoustic  model  simulation.  The  vertical  axis  is  acoustic  transmission  loss  in 
decibels  (referenced  to  1  m).  The  horizontal  axis  is  the  range  in  km.  The  acoustic  source  is  at  a  frequency  of 
400  Hz,  and  the  receiver  depth  is  90  m.  The  transmission  loss  for  the  center  track  is  represented  by  the  gray 
curve.  The  dotted  line  is  transmission  loss  for  the  parallel  track  that  is  4  km  to  the  west.  The  solid  black  line  is 
the  parallel  track  4  km  to  the  east.  The  three  transmission  loss  curves  show  differences  over  the  whole  range 
with  the  most  notable  at  about  40  km  where  differences  of  nearly  20  dB  arc  evident.  This  large  discrepancy  is 
due  to  the  different  bathymetry  along  each  track.  The  strait  is  located  around  40  km  in  range  and  the  bathym¬ 
etry  depth  at  that  range  is  less  for  the  East  and  West  tracks  than  for  the  Center  track.  Thus,  to  a  large  degree, 
the  differences  in  the  bathymetry  are  contributing  to  the  differences  in  transmission  loss  between  the  three 
tracks. 

To  isolate  the  effects  of  the  3D  sound  speed  field,  the  bathymetry  was  made  the  same  for  all  three  tracks. 
Center.  West  and  East.  The  Center  track  bathymetry  was  used  for  all  three  tracks.  The  results  of  this  set  of 
simulations  is  showr  in  Fig.  17.  With  the  bathymetry  the  same  for  all  three  tracks,  the  only  environmental  dif¬ 
ference  left  is  due  the  three  dimension  effects  of  the  sound  speed  field.  The  transmission  loss  in  Fig.  17  shows 
that  the  environments  tor  the  three  tracks  are  acoustically  similar  up  to  a  range  of  about  14  km,  where  the 
effects  of  the  difference  in  the  sounds  speed  fields  are  beginning  to  modify  the  transmission  loss.  After  the 
acoustic  track  passes  the  sill  the  differences  become  even  more  pronounced. 

This  simple  acoustic  simulation  helps  to  demonstrate  the  acoustic  impact  of  accounting  for  the  3D  aspects 
of  the  oceanography.  Ignoring  the  variation  in  the  sound  speed  field  can  produce  significant  differences  in  the 
acoustic  propagation. 


7*  Conclusions 

Three-dimensional  simulations  of  the  Strait  of  Messina  and  Gulf  of  Gioia  region  were  conducted  with  the 
newly  developed  and  documented  ocean  version  of  non  hydrostatic  model  EULAG,  proven  successful  in  sim¬ 
ulating  a  variety  of  rotating  stratified  flows.  The  simulation  were  initialized  from  measured  temperature  and 
salinity  profiles  and  forced  by  the  existing  semidiurnal  tidal  magnitudes  obtained  from  a  tidal  model.  The  sim¬ 
ulation  showed  that  as  the  semidiurnal  tide  moves  over  the  sill  in  the  Strait  of  Messina,  it  generates  a  depres¬ 
sion.  The  depression  results  in  right  and  left  propagating  waves  that  form  internal  bores.  During  tidal 
reversals,  the  formed  bores  undergo  hydraulic  jumps  over  the  sill  and  salinity  and  temperature  fronts  develop. 
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The  fronts  are  defined  by  the  internal  bore  boundaries.  Over  the  sill  the  interface  formed  by  the  front  separates 
water  masses.  The  interface  can  be  located  anywhere  between  the  surface  and  the  sill. 

The  internal  bores  that  propagate  away  from  the  sill,  steepen  on  the  leading  edge  through  nonlinear  effects 
and  disintegrate  into  solitary  wave  trains  as  amplitude  and  frequency  dispersion  set  in.  This  has  been  shown 
previously  with  a  2D  layer  model  by  Brand  et  al.  (1997).  The  propagating  solitary  wave  trains  towards  the 
Gulf  of  Gioia  exhibit  curvature  effects  in  horizontal  planes.  In  the  Gulf  of  Gioia,  the  solitary  wave  trains  shoal. 
In  the  shoaling  process  a  depression  is  formed  along  the  Gulf  of  Gioia  slope.  The  depression  pushes  less  dense 
surface  water  to  depths  of  about  150  m.  Moorings  placed  in  the  model  predictions  indicate  the  presence  of 
lighter  water  at  depths  below  the  surface  during  the  shoaling  of  solitary  wave  trains  in  the  Gulf  of  Gioia.  This 
suggests  mixing  of  water  masses  at  depth.  Mixed  patches  of  water  have  been  observed  by  Guardiani  et  al. 
(1988),  in  the  Gulf  of  Gioia,  during  the  JANE  1984  cruise. 

The  model  results  and  data  were  compared  by  conducting  a  wavelet  analysis.  The  wavelengths  were  tracked 
by  the  spines  (maximum  intensity  for  each  wavelength)  at  various  times.  From  the  slope  of  the  variations, 
phase  speed  were  derived  as  a  function  of  wavelength.  For  the  parameters  extracted  from  CTD  measurements 
and  existing  tidal  conditions,  phase  speed  distributions  for  wavelengths  ranging  from  about  0.6  m  to  1.6  km 
were  obtained.  A  correction  for  feature  distortion  was  applied  to  the  last  tow  of  the  measured  soliton  phase 
speed  distribution.  The  model  predicted  phase  speed  versus  wavelength  had  similar  trends  to  the  phase  speed 
derived  from  data.  The  model  predicted  phase  speed  magnitudes  ranged  from  0.88  m  s  1  to  0.93  ms1.  The 
phase  speeds  derived  from  data  ranged  from  0.77  ms  1  to  0.88  ms  Model  predicted  wavelengths  and 
amplitudes  were  within  15%  of  CTD  chain  measurements. 

Parameters  of  thermocline  and  halocline  depths  were  varied  from  the  measured  conditions  at  the  CTD 
station  by  increasing  them.  The  semidiurnal  tidal  forcing  was  also  varied  by  reducing  it  to  half  of  the  extracted 
barotropic  tidal  forcing  at  the  boundary.  A  variation  of  the  tidal  parameter  to  half,  yielded  a  lower  model 
predicted  phase  velocity  distribution  versus  wavelength.  The  magnitude  of  the  phase  velocity  ranged  from 
about  0.805  ms1  to  0.75  m  s_!  and  exhibited  a  different  wavelength  dependence  from  the  larger  tidal  param¬ 
eter  case.  An  increase  of  the  thermocline  and  halocline  depth  parameters  yielded  a  higher  phase  speed  versus 
wavelength  for  wavelengths  greater  than  1  km.  The  phase  speed  as  a  function  of  wavelength  had  an  increasing 
trend  as  the  wavelength  increased. 

Finally,  the  acoustic  predictions  using  a  3D  ocean  environment  can  be  significantly  different  when  com¬ 
pared  with  similar  acoustic  predictions  using  a  2D  ocean  environment.  In  our  simulations,  even  with  bathym¬ 
etry  effects  removed,  5  dB  differences  in  acoustic  signals  were  observed  over  a  range  of  70  km. 
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Appendix.  Numerical  apparatus 

Given  (7),  the  prognostic  Eqs.  ( 3)— (5)  can  be  written  in  the  symbolic  form  of  an  Eulerian  conservation  law 

^+V(p-v»  =  p-/?,  (16) 

where  \j/  symbolizes  components  of  v,  p\  or  s  and  R  denotes  the  associated  rhs.  The  archetype  problem  (16)  is 
approximated  to  second-order  accuracy  in  space  and  time  using  an  NFT  algorithm 

^+1  =^,(iA)  +  0.5Jf^+1  =ij/i  +  0.5AtRHl+l;  (17) 

where  ^"+l  is  the  solution  sought  at  the  grid  point  (?+l,X|),  ij/  =  ij/"  +  0.5 AtR".  and  sst  denotes  a  two-time-lcvel 
flux-form  Eulerian  nonoscillatory  advective-transport  operator  MPDATA;  sec  Smolarkicwicz  (2006)  for  a 
recent  overview. 
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Eq.  (17)  represents  a  system  implicit  with  respect  to  all  dependent  variables  in  (3)  ami  (4),  because  all  forc¬ 
ing  terms  are  assumed  to  be  unknown  at  n  +  1.  Notably,  the  salinity  Eq.  (5)  has  zero  right-hand  side,  where¬ 
upon  s  has  a  status  of  a  passive  scalar.  For  the  physical  velocity  vector  v,  (17)  can  be  written  compactly  as 

Vi  =  Vj  -0,5df{G(V7t')).  +  0,5zhR|(v,  p*),  (18) 

where 

Ri M)  0.5J/((GTv)  ■  VpJ.)  -  (f  x  (v  -  vc)),  (19) 

P  0 

accounts  for  the  implicit  representation  of  the  buoyancy  via  (4),  and  the  superscript  n  4-  1  has  been  dropped  as 
there  is  no  ambiguity.  Given  a  grid  collocated  with  respect  to  all  prognostic  variables,  (18)  can  be  inverted 
locally  to  construct  expressions  for  the  solenoidal  velocity  components  that  are  subsequently  substituted  into 
(2)  to  produce  an  elliptic  equation  for  pressure 

« 0,  (20) 

where  G 1  v  —  (I  -  0,5zl/R)_lG(  W)J  =  vs  defined  in  (6),  and  the  additional  hats  and  prime  denote  straight¬ 
forward  algebraic  modifications;  cf.  Prusa  and  Smolarkiewicz  (2003)  for  the  complete  exposition.  Boundary 
conditions  imposed  on  v*  n,  subject  to  the  integrability  condition  j^p'vs  ■  ndcr  =  0,  imply  the  appropriate 
boundary  conditions  on  it"  (Prusa  and  Smolarkiewicz,  2003;  Wedi  and  Smolarkicwicz,  2004).  The  resulting 
boundary  value  problem  is  solved  using  a  preconditioned  generalized  conjugate  residual  GCR(A)  algorithm 
(Eisenstat  el  al.,  1983);  for  further  discussion,  see  Smolarkiewicz  et  ah  (2004)  and  references  therein.  Given 

the  updated  pressure,  and  hence  the  updated  solenoidal  velocity,  the  updated  physical  and  contra  variant 

velocity  components  are  constructed  from  the  solenoidal  velocities  using  transformations  (8)  and  (6), 
respectively. 
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